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Abstract 

We are interested in the finite element solution of elliptic problems 
with a right-hand side of the single layer distribution type. Such prob- 
lems arise when one aims at accounting for a physical hypersurface (or 
line, for bi-dimensional problem), but also in the context of fictitious 
domain methods, when one aims at accounting for the presence of an 
inclusion in a domain (in that case the support of the distribution is 
the boundary of the inclusion). The most popular way to handle nu- 
merically the single layer distribution in the finite element context is to 
spread it out by a regularization technique. An alternative approach 
consists in approximating the single layer distribution by a combina- 
tion of Dirac masses. As the Dirac mass in the right hand side does 
not make sense at the continuous level, this approach raises partic- 
ular issues. The object of the present paper is to give a theoretical 
background to this approach. We present a rigorous numerical analy- 
sis of this approximation, and we present two examples of application 
of the main result of this paper. The first one is a Poisson problem 
with a single layer distribution as a right-hand side and the second one 
is another Poisson problem where the single layer distribution is the 
Lagrange multiplier used to enforce a Dirichlet boundary condition on 
the boundary of an inclusion in the domain. Theoretical analysis is 
supplemented by numerical experiments in the last section. 

1 Introduction 



We are interested in the numerical handling of elliptic problems with a right- 
hand side of the single layer distribution type, i.e. problems of the form 
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—An = (p5y in J7, 
u = on <9f2 



where ft is a bounded domain in M. d , and <^<5 7 , in H 1 (fl), can be formally 
written as 




where 7 is a d — 1 dimensional smooth manifold (typically the boundary of 
a connected subdomain wCCd) and ^ is in a Sobolev space H~ l l 2+S {j) 
for s > 0. 

This type of problem is commonly encountered in numerical modeling. 
Peskin and Mc Queen [2] modeled the motion of the heart wall by a collection 
of elastic fiber immersed in a fluid. The force applied to the fluid by these 
fibers is a single layer distribution. In a similar context, Gerbeau et al. [8] 
modeled more recently a valve by an elastic interface immersed in a fluid. 
Likewise, Cottet and Maitre in [T5] track the interface of an elastic membrane 
immersed in a fluid by a level set method and the elastic forces on the 
membrane appears in the right-hand side of the Navier-Stokes equations. 
In [19] the interface between two different fluids is tracked with a level set 
method and the single layer distribution is the surface tension. 

This kind of problem also appears in the context of fictitious domain 
methods. For examples in |12] and |13j the authors use Lagrange multipliers 
to enforce Dirichlet boundary conditions. These Lagrange multipliers act as 
singular distributions of forces supported by the boundary. In [llj, such 
an approach is proposed to simulate the sedimentation of rigid particles in 
a fluid. The rigid body constraint is enforced with Lagrange multipliers 
supported on the boundary of the rigid bodies. The single layer distribution 
also appears in the so called Fat Boundary Method (see [2], or [7] for 
a full analysis of the method). The resolution of a Poisson like problem 
in a domain with holes with this method consists of splitting the Poisson 
problem into two new problems: a global problem, which is solved on a 
regular mesh, and a local problem around the holes, which introduces the 
single layer distribution on the boundary of the holes. 

There exist several approaches to solve these problems numerically, the 
choice of one of them depending on the general framework of the discrete 
problem. In [2], the authors have chosen to solve the problem in a finite 
difference framework. They take a lattice on the whole domain to solve the 
fluid part of the problem and a uniform collection of points on the fibers, 
which are not necessarily the same as the points of the lattice. In order 
to take into account the force of the fibers to the fluid, they regularize the 
Dirac Delta function on every points on the fibers to spread the single layer 
distribution out on the computational lattice. This is the so called Immersed 
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Boundary (IB) method. Applications of the IB method can be found in 
section 9 of [JJ and, for instance, in [3j[ll[5j[6] for more recent ones. Another 
way to deal with the Dirac Delta functions is to directly inject the jump of 
the normal derivatives of the corresponding solution in the finite difference 
scheme. This method, called the Immersed Interface Method (IIM), has 
been introduced by Leveque and Li in [9] . For applications of the IIM, one 
can look in [TU] and the references therein. 

In a variational context, one can compute the integrals involving the 
single layer distribution of the variational formulation in some ways. It 
is the case in a finite element context as in [12] where the authors use 
a constant piece-wise element to discretize the Lagrange multiplier space 
and then compute the integrals exactly. The IB method has been adapted 
to the finite element framework where one can deal with the Dirac Delta 
function in a variational way (see [16} [T7]). In the same spirit, the idea is 
to discretize the single layer distribution by using a combination of Dirac 
masses at a collection of points describing the interface. That is the method 
used in [8j [Til 03] to approximate the Lagrange multiplier defined on an 
interface of the computational domain. This method requires no additional 
mesh but only the collection of points on the interface and is convenient to 
implement in a code. Moreover the integrals which appear in the variational 
formulation are fast and easy to compute. 

Up to our knowledge, the latter approach has not been justified from 
a theoretical point of view. One of the difficulties is that a well defined 
linear functional, the initial right-hand side which we will consider in iT -1 , 
is replaced by a combination of Dirac masses, which do not make sense at 
the continuous level (as soon as the dimension d is greater than 1). 

The aim of this paper is to give a numerical analysis of this method, 
in the case of a scalar Poisson problem with a single layer distribution as 
a right hand side in the two dimensional setting. This paper is structured 
as follow: In section [2] we establish an error estimate when tp is in H~ l l 2+S 
with < s < 1/2 and show that there is a saturation of the order once s is 
greater than 0.5 (when <p is at least in L 2 ). Then, in section [3] we present 
some examples on which this result can be used to estimate the error. In 
section HI we present the numerical results to validate the error analysis done 
in section [2j 

2 General theorems 

Let Q be a domain in M 2 and let a; be a smooth subdomain with boundary 
7. We consider a single layer distribution supported by 7: 

<p5^ G : v G V = i?o(n) ' — ► (<p, u) H -i/ 2+s(7))Jf i/2- S ( 7 ) , 
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with ip G i? _1 / 2+s (7), < s < 1. The expression above makes sense as 
v G H^(Q), so that its trace on 7 is in H l / 2 (^f) C H l / 2 ~ s . 

Given a conforming triangulation 7X of we denote by Vh the associated 
P 1 finite element space: 

V h = {v h G C°(Tl) , v h]K is affine , MK G T h ) . 

We are interested in approximating t^5 7 over Vj, by a combination of 
Dirac masses. Note that it does not make sense at the continuous level, as 
V = contains no Dirac mass. Approximation properties of such an 

appropriate combination will be expressed by Proposition [TJ It is mainly 
based on the following two lemmas. 

Lemma 1. Let I be the unit interval (0, 1), h > 0, and (Th) a family of 
quasi-uniform triangulations of I. We shall represent Th by its subintervals 
71, . . . , 7n (we drop the explicit dependence 0/7, upon h). Quasi-uniformity 
expresses 

ch < \ < Ch , with < c < C. 

We now consider x%, . . . , xn, with Xi G 7$. For any v G H (I), we denote 
by Vh the corresponding piecewise constant interpolant (see Fig. [7]j 

N 



Vh = y~]v{xj)i~ 



i=l 

For any r, < r < 1/2, one has 

\\v — Vh\\H r < Ch l ~ r \v\ 1 

where \v\ 1 is the H 1 seminorm, and the fractional derivative Sobolev norm 
is defined by 

Ml— l/H'l +1/ / '^-Tg'T . °<,<i/2. 

\y-x\ J 

Proof. Let us start with the integral over I x I. Setting Wh = v — Vh, one 
has 

w h (y) - w h (x)\ 2 _ sr^sr^ f f \v(y) - v(xj) - v{x) + v(xi)\ 2 




EE 



I |l+2r / / I |l+2r 

/•/i ijhihi \y-x\ 



j j \aM ~ t '( x )l 2 j ^T^T / / kfe(y) -w fe (x)| 2 
j J~i%J~i% \y — x \ i ■^ i J~i l J~ij \y — x\ 

A + B. 
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Figure 1: P° interpolant 



The first term (diagonal terms) writes 

\v(y) - v(x)\ 2 



^ = E 



■i \y-x 



l+2r 



< 



\y - X\ 



< 



1 



h 2-2r I I „.M2 



(l-2r) 

As for extradiagonal terms, let us first note that, on any 7$, 

2 



\Wh\' 



r v '{t)dt < h f \v'{t)\ 2 . 



The second term B can then be estimated 



B < 4 



EE/ / 



|l+2r 



j^i ^ J Ti 1^ ~ X l 



2r ' 



j n ->n \y ■ 

In the case r > 0, the last sum above (over j / i) is less than twice 
dy 



h r ' 

dx 







so that 



h \y - x 



l+2r 



< 



2r(l - 2r) 



h 



l-2r 



£?< ^— : ^ 2 " 2r / \v'(t)\ Z . 



r(l - 2r) 



Thus, one obtains 



f f \w h (y) - W h {x)\ 2 C ,2-2r U ,|2 

JlJl ' -< 1+2r lV ^' 



\y-x\ 1+2r ~ r(l-2r) 
The 0-th order term is simply 



/2 ^ — ^ / 1 2 2 i 1 2 



'7. 

Hence the order 1 — r of the lemma. □ □ 



Remarks 

• If r = 0, the .ff r norm is only the I? norm and the result comes from 
the 0-th order term of the proof. 

• This result may be ruled out in the case — 1/2 < r < 0. Indeed, even 
though the semi-norm part of the definition of the H r norm can be 
extended for r < 0, it is not the case for the L 2 norm. One cannot 
get rid of this L 2 norm either, because the constant functions must be 
controlled in a way. Moreover, the integral over / x / is of order h 2 ~ 2r 
even if r < 0, so it must also be true for the expression that replaces 
the 1? norm and controls the constant functions. For instance, taking 
the average of the function is not sufficient because it is still a term of 
order 1. That is why it may be that the order cannot go further than 
1, even for r < 0. 

The second lemma is a trace theorem for discrete functions. 

Lemma 2. Let fl be the unit square in M 2 , and uj CC fl a smooth subdo- 
main, with boundary 7. Let {Th)h be a regular family of triangulations, and 
Vh the associated P 1 Finite Element space. The trace operator 70 maps Vh 
(equipped with H 1 norm) onto H l {~i), with 

C 

l7o(v/i)li )7 < ^1/2 Mi,n • 

Proof. Because Vh is a P 1 function, its gradient is constant in every mesh 
cell of Tfr. Let Q 7 be the set of all the mesh cells which have a non empty 
intersection with 7, one has 

ho(v h )\l,= [\vv h \ 2 = \^ v h\lo, Q \Q^i\ 

J ^ QeQ-y 

= K\l Q \Qnn,\ . 
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Assuming that the radius of curvature of 7 is greater than h, the following 
inequality holds, 

|Qn 7 | < Ch, 



where C is a constant independent of h. Thus, the previous equality writes, 



h 

Qe<2 7 
<r C \ I 2 



which is the estimation of the lemma. 



□ 



□ 



The next proposition is the approximation error of a single layer distri- 
bution by a combination of Dirac masses. 

Proposition 1. Let $7 be the unit square in M 2 , and lo CC $7 a smooth 
subdomain, with boundary 7. Let (T^h be a regular family of triangulations, 
and Vh the associated P 1 Finite Element space. Let < s < 1/2 and 
if G F- 1 / 2+s ( 7 ). 

Let (Sr) be a family of quasi-uniform triangulations of 7. We shall 
represent by its subintervals 71, . . . , 7N- h - 

We now consider x\, . . . , xjy h , with Xi G 7« and the following approxi- 
mation of if, 



1=1 

where the real numbers Xi are chosen as follows, 

Xi = <¥>,l 7i ) • 
There exists a constant C > suc/i i/iat : 



(99, u h ) - <C \l-h s ||^||_ 1/2+Si7 |vft|i 



(1) 



(2) 



(3) 



for all function in and where ( . , . ) denotes the dual pairing between 
F- 1 / 2+s ( 7 ) and H 1 ' 2 -*^). 

Proof. Thanks to the definition of </?^, one has, 



N ~ h 



= {<P,Vh) ~ ( (p,^2v h (xi)l^ 
\ i=l 

< IM|-l/2+s, 7 IK - ^2v h (Xi)l^\\ 1/2 - Sil ■ 



i=l 



7 



The function vr defined by 



v h 



y^v h (xi)i- 



is the piecewise constant interpolant of on 7. Now, in order to use lemma 
[1] we perform a change of variable to map the curve 7 onto [0, 1] , 



X: [0,1] 
t 



7 

X(t) 



The variable t here is the arc length so X{t) is the point on 7 such that 
the length of the arc X(0)X(t) is equal to I7I t. Thus, denoting by w\ the 
function Vh — v% , one has, 



2 

l/2- Sl7 



w h(y) - w h( x ) 



7 j i 

7 1 

Jo 



\y - x\ 



i2-2s 



wloX{t)-w\oX(T) 



i2-2s 



-X'(t)X'(r)dtdT. 



\X{t)-X{r) 

For any £ in [0, 1], is equal to | —y | so the change of variables writes, 

h/l 



l/2-s, 7 



w h h oX 



l/2-s,[0,l] 



The same goes for the L 2 norm. Now, the length of the subintervals 7^ of 
7 is h, so the length of the corresponding subintervals in [0, 1] is h/ |7|. We 
apply the lemma[T]with r = 1/2 — s and h = h/ |7|, 

lk-^lll/2- S ,7 = ItI ll^°^-^ oX lll/2-s,[0,l] 



< 



\i\c 



- | 7 | 

Again, with the change of variables, one has, 

C 



h^\v h oX\ im] 



l^-^lll/2- S ,7 " 



hi 



1/2+s 



~ h l/2+s 



\ v h\ 



1,7 



Now, the trace lemma [2] allows us to extend the H 1 semi-norm on 7 of 
Vh over the whole domain Q, 

C h l / 2 + s . 
IK - ^Hi/2- S)7 < 1/2+fl -f^jr Mi,n , 



which ends the proof. 



□ 



□ 
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Remark 1 : If 1/2 < s < 1, one cannot expect to have a better order 
than the one with s = 1/2. Indeed, here is a one dimensional example where 
(p is a smooth function and the error is still of order 1. 

The domain is (0, 1) and (p is the constant function 1. We take h = h/2 
and the function v h is defined by 



Vh(x) 



h/2 



h 



I 



+ 1 in [0, h] 

elsewhere. 



If we take the point h/2 to approximate (p in [0, h], the error writes, 

{<P, v h) 



<ft v h 




h 
2 



Thus the error is of order one and there is a saturation of the order once s 
is greater than 0.5. 



Remark 2 : In the case s = the definition of the approximation of the 
single layer distribution is different because the characteristic functions are 
not in H 1 ! 2 . Nevertheless, these characteristic functions can be replaced by 
a partition of unity and the result of Proposition Q] holds. Thus, there is still 
convergence if h is such that h/h goes to as h goes to 0. 



3 Applications 

3.1 Error estimate for the Poisson problem 

We consider here a Poisson problem with homogeneous Dirichlet boundary 
conditions and a single layer distribution as a right-hand side. We are inter- 
ested in the error estimate when the single layer distribution is approximated 
by a combination of Dirac masses. 

Proposition 2. Let Q be the unit square and wCCSla smooth subdomain, 
with boundary 7. We consider the following Poisson problem, 

—Ait = (p5ry in Q 
u = on d£l 

where ip is in H^ 1 /^ 13 ^), < s < 1/2. 

Let (T/j)/j be a regular family of triangulations and Vh the associated 
P 1 Finite Element space. We approximate ip by a function according 
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to Eqs. HP and Denoting by Uh the finite element solution of u and 
assuming that h is of the order of h, the error estimate is 



\\u - U h \\i < C \yh\u\2 + ^ S |MI-l/2+s, 7 J , 

where \.\% denotes the H 2 semi-norm. 

Proof. Thanks to Strang's lemma (see [E]), the error estimate is, 



\u — Uh\ <C inf [ it — Vh\ + sup 



Vh€V h w h €V h \\UJh\ 

The first part is the usual error and is of order 1/2 in our case because 
is a non conformal mesh and Vh is the Finite Element space of P 1 function. 
The second part is given by Proposition [TJ Thus, the error estimate writes, 

" - < C ( Vh\u\ 2 + h B \\<p\\-x/2+sn 



In particular, if h is equal to h, the error estimate becomes, 

\\U - U h \\i < C (Vh\ll\2 + /l S |M|-l/2+s, 7 ) , 

which ends the proof. □ □ 



3.2 Error estimate for the approximation of the solution of 
a saddle-point problem 

We consider a Poisson equation in a domain with a hole and homogeneous 
Dirichlet boundary conditions on each boundary. The homogeneous condi- 
tion on the boundary of the hole is enforced by Lagrange multipliers defined 
on the boundary. These Lagrange multipliers are approximated by a com- 
bination of Dirac masses. First we prove an abstract theorem on the error 
approximation when the discrete Lagrange multipliers space is not included 
in the continuous Lagrange multipliers space. Second, we use this theorem 
and Proposition [T] to give an error estimate of the Poisson problem. 

Proposition 3. Let V be a Hilbert space, a be a bounded elliptic bilinear 
form in V and f be in V' . Let A be another Hilbert space, B £ C(V, A) and 
K the kernel of B. We consider the problem of finding u in K such that 
a(u, v) = (/, v) for all v in K , and its sadle-point formulation 

a{u,v) + (B*\,v) = (f,v) \/veV 
O, Bu) =0 V/16 A. 
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Let Vh be a finite dimensional subspace of V and A h be a finite dimensional 
space, not necessarily included in A. Let B~^ be a bounded linear application 
from V to A h , we denote by the approximation of K 

K\ = {v h G V h , ((Btyvh, H ) = V Mi € A^} . 

The approximate saddle-point problem is 

a{u h ,v h ) + ({B*)*\- h ,v h ) = (f,v h ) Vv h eV h 
(^ B t u h) =0 V^GA^. 

We have the following error estimate 

\\u-u h \\ v < (l + — ) i n f \\ Wh - u \\ + L in f \\£ - (b£\* fji-J 
V a J Wh eKb a n- h eA- h V V n V 

where a is the coercivity constant of a and £ is the linear form in V defined 
by 

a(u,v) + (Z,v) = {f,v) VveV. 

Proof. Uh is the solution of the approximate saddle-point problem, so we 
have, 

a{u h ,v h ) = (/, v h ) \/v h G K% . 
For all Wh G K£, we write Vh = Uh — Wh G K^. Thus, 

a(v h ,v h ) = a(u h -w h ,v h ) 

= (f,v h )-a(w h ,v h ). 

As u is the solution of our problem we have, 

a(u,v h ) + (£,v h ) = (f,v h ) . 

Thus, for any fi h in A^, one has, 

a{v h ,v h ) = a(u,v h ) + (£, v h ) - a(w h ,v h ) - {B~v h , n~ h ) 
= a{u-w h ,v h ) + fJ> h ,Vh} ■ 

Taking the absolute value of the previous equality gives us 

a\\v h \\ 2 < \\a\\\\u - w h \\\\v h \\ + ||£ - (b£) mJv'IKII ■ 
But Vh = Uh — Wh so we have 

a\\u h - w h \\ < \\a\\\\u - w h \\ + U - (#~) n h \\ v > . 



11 



Thus, 



< \\u - w h \\ + \\wh ~ Uh\\ 

< ^i + M^| K _ n || + I|| C _^)%-| k , y Wh£K ^ H£ 

which concludes the proof. □ □ 

Now we define the problem we are interested in. Let be a bounded 
domain in M 2 and O a non empty open subset of Q. We denote by V the 
following problem 



(V) 



-Au = f in ft \ O 
u = on 7 
u = on Oft , 



where / is a function of L 2 (S7 \ O) and 7 denotes the boundary of O. We 
use a fictitious domain method and extend the function / by in O. We 
still denote by / this extension and we denote by V the space Hq(Q) and 
by K the constrained space, 

K={v€V; ^ = 0} . 

The Lagrange multipliers space is A = Lq(j), which is the space of L 2 
functions on 7 with zero mean value. Let B be the following application, 

B : V —> A 

v — > v h . 

We have K = ker B and the corresponding saddle-point problem is: 
Find (u, A) G V x A, such that 

Vu-Vv+ Xv = fv \/v G V 
fiu = V/i G A . 



Let ft. > 0, we discretize the boundary 7 by taking points x i ~ h G 7 so that 

the distance between two consecutive points is h. Let C V be a finite 
element space approximation of V. The Lagrange multipliers space and the 
constrained space are approximated as follows: 



A~ h = R N h 



K\ = [v h G V h ; v h (x i:h ) = 0} 
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Here the space is not included in A. The application B becomes 
B\: V h — ► Ar 



"h- Vh ^ 

v h — »■ \Vh{x Lh 



i = l...Nr 

h 



We denote by (V h h ) the following approximated problem: 
Find (v,h, At) G Vh X At, such that 

a{uh,v h ) + (^X H , B^v h \ = (f,v h ) Vv h eV h 

Where a is the bilinear form defined in V x V by 

a(u, v) = Vu • Vw . 
./n 

We can apply the result of Proposition [3] so we have the following error 
estimate, 



\u - Uh\\v < 



1 + inf \\w h -u\\ + - inf M - h*)* fji K \\ . 



As in the previous example, the first part of this error estimate is the usual 
one and the second part is the error done by approximating the linear form 
£ with a combination of Dirac masses. Indeed, for any /it = (/jti)i ( .„ f jvw i n 
A^ and for any x>h in Vh, one has, 



i I \i=l I 

The error estimate of Proposition [3] becomes, 
u — Uh i < 1 H mi wjj -wM — mi 



a / a /^eA,- \\v h \ 

Thus, by applying the result of Proposition [H 



\u-U h \\l<C\ ^H! + ^ ^ 1/2 ||£lk 



Again, if h is equal to h, the error estimate writes, 

\\u - u h \\i < cv^CHi + ||e|| 0,7 ) ■ 
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4 Numerical results 



In this section we present two examples to validate the estimate ([3]) of Propo- 
sition [TJ The first one is a one dimensional example to show that in the case 
h = h the error order is the one predicted by Proposition [TJ The second 
example is a two dimensional numerical example where the solution of a 
Poisson equation with a right-hand side of the single layer distribution type 
is computed. This right-hand side is the Laplacian of a function with a 
jump of its normal derivatives across an interface. The exact solution is 
thus known and we compute the error in the case h = h for various values 
of s. 



4.1 Id example 

Let ip be the function defined by 



1 



nl — S 



where s is positive, and Vh be the affine function defined in [0, h] by 

h 



h 



+ 1. 



The estimate of Proposition [JJ can be computed explicitly. Indeed, one has, 

1 1 1 \ h s+1 



(<P, Vh) 



L 



h y 
(pv h = h 



X s + 1 ) s(s + 1)2 S 2s J h 
Let Xi be in [0, h], the second part is 



<Ph> v h 



¥ V h {Xi) 



™ . h 

A ' 2 



h 



■ h s 
+ 1 I — 

s 



Let < a < 1, we write Xi = ah. Thus, the error estimate writes 

h s+1 



(if, v h ) - (ifi v h ) =C{s) 



h 



where C(s) is a constant depending only on s. Now, if h = Ch, where C is a 
constant, the error is of order s which is the order predicted by Proposition 

m 
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4.2 Two-dimensional example 

In this section f2 is the unit square and 7 is a circle of radius R. To investigate 
the behavior of the approach for a right-hand side ip5^, with <p in 
we build the function tp as an infinite series of sine functions 

00 

ip = C n ~ s sin(nO) 
n=l 

We now build a function u such that the jump of its normal derivatives on 
7 is <p so that u is the solution of a Poisson problem of the form 

-Au = c^<5 7 + / in Q 
u = on c?0 

where / is to be defined later on. 

First we need a cut function to make sure that the solution u is equal to 
zero on the boundary of Let us denote by \e t ne C 2 function defined in 

R by 

( 1 if x < 



XeO) 



-6 . 15 4 10 , 

-^rx b + -rx 4 + 1 if < x < e 







if x > e 



We define p as the distance between 7 and dQ, (see figure [2]). We denote by 
Rmax the sum of R and p. Now, we can define the following function : 



u(r, 9) 



Xp/3 



u [r + ( 1 



R + 

r 

' R 



71=1 



7" \ — n 



i? 



n s 1 sin(rt0) if r > R 



if r < R. 



It is a function which has a jump of its normal derivative on 7 and oscillates 
very fast close to 7. Its Laplacian outside the particle is given by 



Xp/3 



71=1 
OO 



T \ ~ n 



R 



n 



-s-l 



sin(n6 



1 — 2n / r \ -1-1 



71=1 



The jump ip of its normal derivatives is : 

Rm.r, 



1 



R 



R \R 



i? 2 



- 2^ n sm(nfc') 



71=1 



_1 sin(n#) . 



15 



p R Rmax 



x = i 



= 

9. 



Figure 2: Notations 

which is in H~ l l 2+S ~ e {^f) for all e > 0. This means that we have the exact 
solution of the following problem 

— Au = ip5ry + f in £1 
u = on dO, 

where / is given by the Laplacian of u outside the particle. 

To compute the numerical solution of this problem we truncate the series 
at some order N greater than 2ir/h so that there is at least one period of 
the function sin(iV0) in a mesh cell. We compute the numerical solution for 
several s using Q\ finite elements and compare it to the exact solution (see 
figure [3|). In all the tests, the order N is equal to 2 12 and h takes its values 
from 2 -6 to 2 -10 . The space step h is equal to h so that the H 1 error should 
be of the form, 

\\u - u h \\i < C (Vh\u\ 2 + h s \\tp\\_ l / 2+ s,^ ■ 

We can see in figure [3] that the error order seems to stay close to 0.6 
once s is greater than 0.5. This is because of the space approximation of 
order 1/2 since we use a cartesian mesh with finite elements of order 1, but 
it could also be because of the saturation of the error order when s is greater 
than 0.5 (see the remark right after Proposition [TJ . 

Figure H] shows the order of the error as s goes to 1. When s is less 
than 0.5, the slope is equal to 1 which is expected. The numerical order is 
always better than the theoretical one, but this kind of super-convergence 
often happens when one compute numerical solutions of PDE. 
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s = 0.125 (slope = 0.2) 1 

8 = 0.25 (slope = 0.32) >( 

s =0.375 (slope = 0.45) — 

s = 0.5 (slope = 0.56) — [=^- 
s =0.625 (slope = 0.62) 

= 0.75 (slope = 0.63) A 
0.875 (slope = 0.63) 

s = 1 (slope = 0.63) — X— 



Figure 3: H 1 error for some value of s 
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Figure 4: Order of the H 1 error with respect to s 
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